## Project: Autocracy UPR project ----------------------------------------------
##        05 Analyses             ----------------------------------------------
## Author:
##   - Chun-Young Park (UGA)
##   - Sang-Hoon Park (UofSC)
##
## Updated: Aug 31 2024
## R version 4.4.0 (2024-04-24)
## Platform: aarch64-apple-darwin20
## Running under: macOS Sonoma 14.6.1

if (!require("devtools")) install.packages("devtools")
if (!require("pacman")) install.packages("pacman")
if (!require("rqog")) remotes::install_github("ropengov/rqog")
pacman::p_load(rqog, rio, tidyverse, dendextend, splitstackshape, psych, dlookr,
               lmtest, aod, texreg, modelsummary, MASS, tidyverse)
if (!require("dendextendRcpp")) devtools::install_github('talgalili/dendextendRcpp')
if (!require("vdemdata")) devtools::install_github("vdeminstitute/vdemdata")
if (!require("ggrepel")) install.packages("ggrepel")
if (!require("patchwork")) install.packages("patchwork")
if (!require("futurevisions")) devtools::install_github("JoeyStanley/futurevisions")
if (!require("wesanderson")) devtools::install_github("karthik/wesanderson")

dlookr::import_google_font("Barlow Semi Condensed")

theme_clean <- function() {
  theme_minimal(base_family = "Barlow Semi Condensed") +
    theme(panel.grid.minor = element_blank(),
          plot.title = element_text(face = "bold", color = "black"),
          plot.subtitle = element_text(face = "bold", color = "black"),
          axis.title = element_text(family = "Barlow Semi Condensed Medium", 
                                    color = "black", size = rel(1.5)),
          axis.text = element_text(color = "black", size = rel(1.3)),
          strip.text = element_text(family = "Barlow Semi Condensed",
                                    face = "bold", size = rel(1), hjust = 0,
                                    color = "black"),
          strip.background = element_rect(fill = "white", color = NA),
          plot.caption = element_text(hjust = 0, color = "black"),
          legend.position = "bottom")
}
showtext::showtext.auto(T)
ggplot2::theme_set(theme_clean())

### Import analysis data -------------------------------------------------------
upr_vdem_hrs <- readRDS("data/data for analysis/upr_vdem_hrs.RDS")

### Dummies for regime dyads ---------------------------------------------------
#### Main: Lexical Index -------------------------------------------------------
upr_vdem_hrs |> mutate(
  dd_lexi = case_when(
    is.na(regime_dyads_lexi) ~ NA_integer_,
    regime_dyads_lexi == "Democracy Dyads" ~ 1L,
    T ~ 0L),
  da_lexi = case_when(
    is.na(regime_dyads_lexi) ~ NA_integer_,
    regime_dyads_lexi == "Demo(REV)-Auto(SUR) Mixed Dyads" ~ 1L,
    T ~ 0L),
  ad_lexi = case_when(
    is.na(regime_dyads_lexi) ~ NA_integer_,
    regime_dyads_lexi == "Auto(REV)-Demo(SUR) Mixed Dyads" ~ 1L,
    T ~ 0L),
  aa_lexi = case_when(
    is.na(regime_dyads_lexi) ~ NA_integer_,
    regime_dyads_lexi == "Autocracy Dyads" ~ 1L,
    T ~ 0L)) -> sample

#### Robustness check ----------------------------------------------------------
##### Polyarchy Index < 0.5 ----------------------------------------------------
sample |> mutate(
  dd_poly = case_when(
    is.na(regime_dyads) ~ NA_integer_,
    regime_dyads == "Democracy Dyads" ~ 1L,
    T ~ 0L),
  da_poly = case_when(
    is.na(regime_dyads) ~ NA_integer_,
    regime_dyads == "Demo(REV)-Auto(SUR) Mixed Dyads" ~ 1L,
    T ~ 0L),
  ad_poly = case_when(
    is.na(regime_dyads) ~ NA_integer_,
    regime_dyads == "Auto(REV)-Demo(SUR) Mixed Dyads" ~ 1L,
    T ~ 0L),
  aa_poly = case_when(
    is.na(regime_dyads) ~ NA_integer_,
    regime_dyads == "Autocracy Dyads" ~ 1L,
    T ~ 0L)) -> sample

##### Boix, Rosato, and Miller -------------------------------------------------
sample |> mutate(
  dd_brm = case_when(
    is.na(regime_dyads_brm) ~ NA_integer_,
    regime_dyads_brm == "Democracy Dyads" ~ 1L,
    T ~ 0L),
  da_brm = case_when(
    is.na(regime_dyads_brm) ~ NA_integer_,
    regime_dyads_brm == "Demo(REV)-Auto(SUR) Mixed Dyads" ~ 1L,
    T ~ 0L),
  ad_brm = case_when(
    is.na(regime_dyads_brm) ~ NA_integer_,
    regime_dyads_brm == "Auto(REV)-Demo(SUR) Mixed Dyads" ~ 1L,
    T ~ 0L),
  aa_brm = case_when(
    is.na(regime_dyads_brm) ~ NA_integer_,
    regime_dyads_brm == "Autocracy Dyads" ~ 1L,
    T ~ 0L)) -> sample

##### PolityIV < 7 -------------------------------------------------------------
sample |> mutate(
  dd_pol4 = case_when(
    is.na(regime_dyads_pol4) ~ NA_integer_,
    regime_dyads_pol4 == "Democracy Dyads" ~ 1L,
    T ~ 0L),
  da_pol4 = case_when(
    is.na(regime_dyads_pol4) ~ NA_integer_,
    regime_dyads_pol4 == "Demo(REV)-Auto(SUR) Mixed Dyads" ~ 1L,
    T ~ 0L),
  ad_pol4 = case_when(
    is.na(regime_dyads_pol4) ~ NA_integer_,
    regime_dyads_pol4 == "Auto(REV)-Demo(SUR) Mixed Dyads" ~ 1L,
    T ~ 0L),
  aa_pol4 = case_when(
    is.na(regime_dyads_pol4) ~ NA_integer_,
    regime_dyads_pol4 == "Autocracy Dyads" ~ 1L,
    T ~ 0L)) -> sample

##### V-Dem RoW  -------------------------------------------------------------
sample |> mutate(
  dd_vdem = case_when(
    is.na(regime_dyads_vdem) ~ NA_integer_,
    regime_dyads_vdem == "Democracy Dyads" ~ 1L,
    T ~ 0L),
  da_vdem = case_when(
    is.na(regime_dyads_vdem) ~ NA_integer_,
    regime_dyads_vdem == "Demo(REV)-Auto(SUR) Mixed Dyads" ~ 1L,
    T ~ 0L),
  ad_vdem = case_when(
    is.na(regime_dyads_vdem) ~ NA_integer_,
    regime_dyads_vdem == "Auto(REV)-Demo(SUR) Mixed Dyads" ~ 1L,
    T ~ 0L),
  aa_vdem = case_when(
    is.na(regime_dyads_vdem) ~ NA_integer_,
    regime_dyads_vdem == "Autocracy Dyads" ~ 1L,
    T ~ 0L)) -> sample

##### Skaaning LIED  -------------------------------------------------------------
##### LIED < 5
sample |> mutate(
  dd_lied5 = case_when(
    is.na(regime_dyads_lexi5) ~ NA_integer_,
    regime_dyads_lexi5 == "Democracy Dyads" ~ 1L,
    T ~ 0L),
  da_lied5 = case_when(
    is.na(regime_dyads_lexi5) ~ NA_integer_,
    regime_dyads_lexi5 == "Demo(REV)-Auto(SUR) Mixed Dyads" ~ 1L,
    T ~ 0L),
  ad_lied5 = case_when(
    is.na(regime_dyads_lexi5) ~ NA_integer_,
    regime_dyads_lexi5 == "Auto(REV)-Demo(SUR) Mixed Dyads" ~ 1L,
    T ~ 0L),
  aa_lied5 = case_when(
    is.na(regime_dyads_lexi5) ~ NA_integer_,
    regime_dyads_lexi5 == "Autocracy Dyads" ~ 1L,
    T ~ 0L)) -> sample

sample |> mutate(
  dd_lied4 = case_when(
    is.na(regime_dyads_lexi4) ~ NA_integer_,
    regime_dyads_lexi4 == "Democracy Dyads" ~ 1L,
    T ~ 0L),
  da_lied4 = case_when(
    is.na(regime_dyads_lexi4) ~ NA_integer_,
    regime_dyads_lexi4 == "Demo(REV)-Auto(SUR) Mixed Dyads" ~ 1L,
    T ~ 0L),
  ad_lied4 = case_when(
    is.na(regime_dyads_lexi4) ~ NA_integer_,
    regime_dyads_lexi4 == "Auto(REV)-Demo(SUR) Mixed Dyads" ~ 1L,
    T ~ 0L),
  aa_lied4 = case_when(
    is.na(regime_dyads_lexi4) ~ NA_integer_,
    regime_dyads_lexi4 == "Autocracy Dyads" ~ 1L,
    T ~ 0L)) -> sample
sample <- sample |> janitor::clean_names()
saveRDS(sample, "data/data for analysis/sample.RDS")
rio::export(sample, "data/data for analysis/sample.dta")

### Model Estimations ------------------------------------------------------------------------------
### Due to the model and data constraints, the estimation processes take a while.
### If you want to check the model outputs, just import the model objects 
### from the data/data for analysis folder

#### Ordered Logistic Models -----------------------------------------------------------------------
##### Model 1: Baseline model ----------------------------------------------------------------------
sample_df <- as.data.frame(sample)

ologit_model1 <- MASS::polr(
  as.factor(severity) ~ 
    aa_lexi + da_lexi + ad_lexi +
    as.factor(to_cow) + as.factor(from_cow) +
    as.factor(year),
  data = sample_df, Hess = T,
  method = "logistic")

gc()

### Some dropped coefficients are from reviewer- or states under review-fixed effects.
texreg::screenreg(ologit_model1, omit.coef = "as.factor")

##### Model 2: Full model --------------------------------------------------------------------------
ologit_model2 <- MASS::polr(
  as.factor(severity) ~ 
    aa_lexi + da_lexi + ad_lexi +
    idealdistance +
    politysenderminustarget + 
    alliance +
    share_combined +
    reviewer_revieweed +
    to_mem + from_mem + regionalism +
    sender_hrs + target_hrs +
    migrantscat + socioecon + vulnerable +
    womenchild + physint + justicecat + political +
    as.factor(to_cow) + as.factor(from_cow) +
    as.factor(year),
  data = sample_df, Hess = T,
  method = "logistic")
gc()
### Some dropped coefficients are from reviewer- or states under review-fixed effects.
texreg::screenreg(ologit_model2, omit.coef = "as.factor")

##### With Bilateral Trade Variables
##### Model 3: Full model --------------------------------------------------------------------------

ologit_model3 <- MASS::polr(
  as.factor(severity) ~ 
    aa_lexi + da_lexi + ad_lexi +
    politysenderminustarget + 
    lagtrade1 +
    alliance + idealdistance + 
    reviewer_revieweed +
    to_mem + from_mem + regionalism +
    sender_hrs + target_hrs +
    migrantscat + socioecon + vulnerable +
    womenchild + physint + justicecat + political +
    as.factor(to_cow) + as.factor(from_cow) +
    as.factor(year),
  data = sample_df, Hess = T,
  method = "logistic")
gc()

ologit_model4 <- MASS::polr(
  as.factor(severity) ~ 
    aa_lexi + da_lexi + ad_lexi +
    politysenderminustarget + 
    lagtrade3 +
    alliance + idealdistance + 
    reviewer_revieweed +
    to_mem + from_mem + regionalism +
    sender_hrs + target_hrs +
    migrantscat + socioecon + vulnerable +
    womenchild + physint + justicecat + political +
    as.factor(to_cow) + as.factor(from_cow) +
    as.factor(year),
  data = sample_df, Hess = T,
  method = "logistic")
gc()

texreg::screenreg(list(ologit_model1, ologit_model2, 
                       ologit_model3, ologit_model4),
                  omit = "as.factor")

saveRDS(ologit_model1, "data/estimation/ologit_model1.RDS");
saveRDS(ologit_model2, "data/estimation/ologit_model2.RDS");
saveRDS(ologit_model3, "data/estimation/ologit_model3.RDS");
saveRDS(ologit_model4, "data/estimation/ologit_model4.RDS")


#### Pull out the results as .tex
texreg::texreg(list( 
  ologit_model1, ologit_model2),
  omit.coef = "as.factor",
  custom.model.names=c(
    "Model 1",
    "Model 2"),
  reorder.coef = c(1, 2, 3, 9, 6, 8, 7, 10, 11, 12, 13, 15, 14,
                   16, 17, 18, 19, 20, 21, 22, 4, 5),
  custom.coef.map = list(
    "aa_lexi" = "Autocracy Dyads",
    "da_lexi" = "Demo$_{\\mathrm{Reviewer}}$-Auto$_{\\mathrm{SuR}}$ Mixed Dyads", 
    "ad_lexi" = "Auto$_{\\mathrm{Reviewer}}$-Demo(SUR)",
    "1|2" = "(Cut 1)",
    "2|3" = "(Cut 2)",
    "politysenderminustarget" = "Joint Democracy", 
    "share_combined" = "Share of Combined GDP",
    "alliance" = "Alliance", 
    "idealdistance" =  "Geopolitical Affinity",
    "reviewer_revieweed" = "Reviewed Reviewer",
    "to_mem" = "HRC Member$_{\\mathrm{SuR}}$", 
    "from_mem" = "HRC Member$_{\\mathrm{Reviewer}}$",
    "regionalism"= "Same Region",
    "sender_hrs" = "Human Rights Score$_{\\mathrm{Reviewer}}$",
    "target_hrs" = "Human Rights Score$_{\\mathrm{SuR}}$",
    "migrantscat"= "Migration", 
    "socioecon" = "Socio-Economic Rights", 
    "vulnerable" = "Vulnerable Populations",
    "womenchild" =   "Women, Children \\& Trafficking", 
    "physint" = "Physical Integrity Rights",
    "justicecat" =   "Justice", 
    "political" = "Speech \\& Political Participation"),
  single.row = TRUE,
  include.rsquared=FALSE, 
  custom.gof.rows=list(`Country-fixed` = c("Yes", "Yes"),
                       `Year-fixed` = c("Yes", "Yes")),
  stars=c(0.001, 0.01, 0.05), 
  #custom.note="%stars. Year effects not significant.",
  caption=
    "\\label{tab1} Ordered logit models on UPR recommendations shaming behaviors",
  file = "tables/table3.tex")

#### Likelihood ratio test (compared to the intercept-only model) --------------
ologit_model_null1 <- 
  polr(ologit_model1$model[, 1] ~ 1, method = "logistic",
       Hess = TRUE)
lrtest(ologit_model1, ologit_model_null1)

ologit_model_null2 <- 
  polr(ologit_model2$model[, 1] ~ 1, method = "logistic",
       Hess = TRUE)
lrtest(ologit_model2, ologit_model_null2)

#### Testing the Cut Points ----------------------------------------------------
##### H0: Cut1 = Cut2
wald.test(Sigma = vcov(ologit_model1), 
          b = c(coef(ologit_model1), ologit_model1$zeta),
          L = t(as.matrix(c(rep(0, 348), 1, -1))))

#### Wald test: -----------------------------------------------------------------
### Chi-squared test:
### X2 = 35508.8, df = 1, P(> X2) = 0.0
### The results show that we have enough empirical evidence to reject 
### the null hypothesis that the two cut points are not different. 
### It implies that the ranks in the dependent variable have distinctive features
### against each other.

#### First Differences in Predicted Probabilities -------------------------------
#### Simulate coefficients and cut points (King, Tomz, and Wittenberg 2000)
set.seed(1234)

simb_ologit1 <- mvrnorm(n = 4000, 
                        mu = c(coef(ologit_model1), ologit_model1$zeta),
                        Sigma = vcov(ologit_model1))
cut_ologit1 <- cbind(-Inf, simb_ologit1[, 349:350], Inf)

simb_ologit2 <- mvrnorm(n = 4000, 
                        mu = c(coef(ologit_model2), ologit_model2$zeta),
                        Sigma = vcov(ologit_model2))
cut_ologit2 <- cbind(-Inf, simb_ologit2[, 360:361], Inf)

simb_ologit3 <- mvrnorm(n = 4000, 
                        mu = c(coef(ologit_model3), ologit_model3$zeta),
                        Sigma = vcov(ologit_model3))
cut_ologit3 <- cbind(-Inf, simb_ologit3[, 332:333], Inf)

simb_ologit4 <- mvrnorm(n = 4000, 
                        mu = c(coef(ologit_model4), ologit_model4$zeta),
                        Sigma = vcov(ologit_model4))
cut_ologit4 <- cbind(-Inf, simb_ologit3[, 332:333], Inf)

#### Create the baseline profile + Calculate the Predicted Probs

x_dd_ologit2 <- cbind(0, 0, 0, 
                      mean(sample$idealdistance, na.rm = T),
                      mean(sample$politysenderminustarget, na.rm = T),
                      0, 
                      mean(sample$share_combined, na.rm = T),
                      0, 0, 0, 0,
                      mean(sample$sender_hrs, na.rm = T),
                      mean(sample$target_hrs, na.rm = T),
                      0, 0, 0, 0, 0, 0, 0)

xddb_ologit2 <- simb_ologit2[, 1:20] %*% t(x_dd_ologit2)

pr0_ologit2 <- matrix(NA, nrow = 4000, ncol = 3)
for (i in 1:3){ # for each value of the dependent variable...
  pr0_ologit2[, i] <- plogis(cut_ologit2[, i + 1] - xddb_ologit2) - plogis(cut_ologit2[, i] - xddb_ologit2)
}

#### Matrix to store the results
out.mat.ologit2 <- matrix(NA, nrow = 3 * 3, ncol = 3)
rownames(out.mat.ologit2) <- paste0(rep(c("Autocratic Dyads",
                                  "Demo(REV)-Auto(SUR) Mixed Dyads",
                                  "Auto(REV)-Demo(SUR) Mixed Dyads"), 
                                each = 3), "-",
                            rep(1:3, times = 3))
colnames(out.mat.ologit2) <- c("Mean", "Lower", "Upper")

#### Then compute the first differences in each covariates
#### DD -> AA
x_aa_ologit2 <- x_dd_ologit2
x_aa_ologit2[1] <- 1

x_aab_ologit2 <- simb_ologit2[, 1:20] %*% t(x_aa_ologit2)
for (i in 1:3){ # for each value of the dependent variable...
  pr1_ologit2 <- plogis(cut_ologit2[, i + 1] - x_aab_ologit2) - plogis(cut_ologit2[, i] - x_aab_ologit2)
  fdiff_ologit2 <- pr1_ologit2 - pr0_ologit2[, i]
  out.mat.ologit2[i, 1] <- mean(fdiff_ologit2)
  out.mat.ologit2[i, 2:3] <- quantile(fdiff_ologit2, probs = c(0.025, 0.975))
}

#### DD -> DA
x_da_ologit2 <- x_dd_ologit2
x_da_ologit2[2] <- 1
xdab_ologit2 <- simb_ologit2[, 1:20] %*% t(x_da_ologit2)
for (i in 1:3){ # for each value of the dependent variable...
  pr1_ologit2 <- plogis(cut_ologit2[, i + 1] - xdab_ologit2) - plogis(cut_ologit2[, i] - xdab_ologit2)
  fdiff_ologit2 <- pr1_ologit2 - pr0_ologit2[, i]
  out.mat.ologit2[i + 3, 1] <- mean(fdiff_ologit2)
  out.mat.ologit2[i + 3, 2:3] <- quantile(fdiff_ologit2, probs = c(0.025, 0.975))
}

#### DD -> AD
x_ad_ologit2 <- x_dd_ologit2
x_ad_ologit2[3] <- 1
xadb_ologit2 <- simb_ologit2[, 1:20] %*% t(x_ad_ologit2)
for (i in 1:3){ # for each value of the dependent variable...
  pr1_ologit2 <- plogis(cut_ologit2[, i + 1] - xadb_ologit2) - plogis(cut_ologit2[, i] - xadb_ologit2)
  fdiff_ologit2 <- pr1_ologit2 - pr0_ologit2[, i]
  out.mat.ologit2[i + 6, 1] <- mean(fdiff_ologit2)
  out.mat.ologit2[i + 6, 2:3] <- quantile(fdiff_ologit2, probs = c(0.025, 0.975))
}
out.mat.ologit2 # print the simulation results

# Visualize the results
library(wesanderson)
out_mat_ologit2 <- as.data.frame(out.mat.ologit2)
out_mat_ologit2 |> mutate(
  sig = if_else(Lower < 0 & Upper > 0, 
                "Insignificant", "Significant")
) -> out_mat_ologit2
out_mat_ologit2$sig <- factor(out_mat_ologit2$sig, 
                              levels = c("Insignificant", "Significant"))
out_mat_ologit2$regimedyads <- rep(1:3, times = 3)

g1_main <- ggplot(out_mat_ologit2[1:3,]) +
  geom_hline(aes(yintercept = 0), linetype = "dashed") +
  geom_pointrange(aes(x = regimedyads, y = Mean, 
                      ymin = Lower, ymax = Upper,
                      color = sig), show.legend = F, size = 1.1) +
  labs(subtitle = "Democratic\nto Autocratic Dyads",
       y = "First Differences in\nPredicted Probabilities",
       x = "\n") + 
  scale_y_continuous(labels = scales::percent_format(),
                     limits = c(-0.25, 0.15)) +
  scale_x_continuous(breaks = c(1:3),
                     labels = c("Praising",
                                "Neutral",
                                "Shaming")) + 
  ggrepel::geom_text_repel(aes(x = regimedyads, y = Mean, 
                               label = paste0(round(Mean, 4)*100, "%"),
                               color = sig),
                           size = 4.5, angle = 90,
                           nudge_x = -0.3,
                           xlim = c(-Inf, Inf), ylim = c(-Inf, Inf),
                           show.legend = F) + 
  scale_color_manual(values = c(wesanderson::wes_palette("Royal1")[2],
                                wesanderson::wes_palette("Royal1")[1])) + 
  theme(legend.title = element_blank(),
        plot.margin = unit(c(0,1,0,0), "cm"))

g2_main <- ggplot(out_mat_ologit2[4:6,]) +
  geom_hline(aes(yintercept = 0), linetype = "dashed") +
  geom_pointrange(aes(x = regimedyads, y = Mean, 
                      ymin = Lower, ymax = Upper,
                      color = sig), show.legend = F, size = 1.1) +
  labs(subtitle = "Democratic\nto Auto(REV)-Demo(SUR) Mixed Dyads",
       y = NULL, x = "\n") + 
  scale_y_continuous(labels = scales::percent_format(),
                     limits = c(-0.25, 0.15)) +
  scale_x_continuous(breaks = c(1:3),
                     labels = c("Praising",
                                "Neutral",
                                "Shaming")) + 
  ggrepel::geom_text_repel(aes(x = regimedyads, y = Mean, 
                               label = paste0(round(Mean, 4)*100, "%"),
                               color = sig),
                           nudge_x = -0.3,
                           size = 4.5, angle = 90,
                           xlim = c(-Inf, Inf), ylim = c(-Inf, Inf),
                           show.legend = F) + 
  scale_color_manual(values = c(wesanderson::wes_palette("Royal1")[1],
                                wesanderson::wes_palette("Royal1")[2])) + 
  theme(legend.title = element_blank(),
        axis.text.y.left = element_blank(),
        plot.margin = unit(c(0,1,0,0), "cm"))

g3_main <- ggplot(out_mat_ologit2[7:9,]) +
  geom_hline(aes(yintercept = 0), linetype = "dashed") +
  geom_pointrange(aes(x = regimedyads, 
                      y = Mean, ymin = Lower, ymax = Upper,
                      color = sig), show.legend = F, size = 1.1) +
  labs(subtitle = "Democratic\nto Demo(REV)-Auto(SUR) Mixed Dyads",
       x = "\n", y = NULL) + 
  scale_y_continuous(labels = scales::percent_format(),
                     limits = c(-0.25, 0.15)) +
  scale_x_continuous(breaks = c(1:3),
                     labels = c("Praising",
                                "Neutral",
                                "Shaming")) + 
  ggrepel::geom_text_repel(aes(x = regimedyads, y = Mean, 
                               label = paste0(round(Mean, 4)*100, "%"),
                               color = sig),
                           size = 4.5, angle = 90,
                           nudge_x = -0.3,
                           xlim = c(-Inf, Inf), ylim = c(-Inf, Inf),
                           show.legend = F) + 
  scale_color_manual(values = c(wesanderson::wes_palette("Royal1")[1],
                                wesanderson::wes_palette("Royal1")[2])) + 
  theme(legend.title = element_blank(),
        axis.text.y.left = element_blank(),
        plot.margin = unit(c(0,1,0,0), "cm"))


g1_main +  g2_main + g3_main + plot_layout(ncol = 3) +
  plot_annotation(
    caption = str_wrap("Note: This figure estimates the change in predicted probability for each value of the dependent variable (praising, neutral, and shaming) as a function of changing the different regime dyads compared to democracy dyads, while holding other covariates at their mean values.
    The y-axis represents the estimated percentage change, with a range from -10% to 10%. 
    Each panel displays the effects for three types of dyads: Autocratic Dyads, Auto(REV)-Demo(SUR) Mixed Dyads, and Demo(REV)-Auto(SUR) Mixed Dyads. 
                       The dashed horizontal line at zero indicates no change, while points above or below represent positive or negative effects, respectively. The grey colors differentiate between the dyads, with shapes indicating statistical insignificance.", 
                       120, exdent = 10)) &
  theme(legend.title = element_blank(),
        strip.text = element_text(size = rel(1.4)),
        axis.title = element_text(size = rel(1.3)),
        axis.text = element_text(size = rel(1.2)),
        legend.text = element_text(size = rel(1.3)),
        legend.margin=margin(0,0,0,0),
        legend.box.margin=margin(0, 0, 20, 0),
        plot.caption = element_text(size = rel(1.4)))

ggsave("figures/fig2.pdf", height = 6, width = 12, dpi = 1600)
